######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate results in Table C.5
######################################################################
rm(list = ls())

#Other packages used
library(tidyverse)
library(lmtest)
library(sandwich)
library(MASS)
library(stargazer)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

###################################
#Produce Table C.5
###################################
#-----------------------
#(1) Load data
#-----------------------
load("data/data_condition.RData")

#-----------------------
#(2) Setup
#-----------------------
#covariates
covs = c("gdppc_log_lag",
         "population_log_lag",
         "election_either_lag",
         "fdi_gdp_lag",
         "debt_gni_lag",
         "oda_gni_lag",
         "polity2_lag",
         "unsc_lag",
         "IdealPointDistance_lag",
         "totalamt_log")

#Treatment
treat.use = "aiib_founder_2016"

#FE
fe.use = c("factor(iso2c)", "factor(year)")

#Data
df.use = df.condition

#-----------------------
#(3) Column (1)
#-----------------------
#DV
dv.use = "pa_count"

#Equation
eq = make.equation(dv.use, treat.use, fe.use)

#Estimate
fit.tab.c5.1 = glm.nb(eq, df.use, control = glm.control(maxit = 20))
se.tab.c5.1 = fit.tab.c5.1 %>%
  coeftest(.,
           vcov = vcovHC(., type = "HC0", cluster = df.use$iso2c))

#-----------------------
#(4) Column (2)
#-----------------------
#DV
dv.use = "sector_count"

#Equation
eq = make.equation(dv.use, treat.use, fe.use)

#Estimate
fit.tab.c5.2 = glm.nb(eq, df.use, control = glm.control(maxit = 20))
se.tab.c5.2 = fit.tab.c5.2 %>%
  coeftest(.,
           vcov = vcovHC(., type = "HC0", cluster = df.use$iso2c))

#-----------------------
#(5) Column (3)
#-----------------------
#DV
dv.use = "pa_count"

#Equation
eq = make.equation(dv.use, treat.use, c(covs, fe.use))

#Estimate
fit.tab.c5.3 = glm.nb(eq, df.use, control = glm.control(maxit = 20))
se.tab.c5.3 = fit.tab.c5.3 %>%
  coeftest(.,
           vcov = vcovHC(., type = "HC0", cluster = df.use$iso2c))

#-----------------------
#(6) Column (4)
#-----------------------
#DV
dv.use = "sector_count"

#Equation
eq = make.equation(dv.use, treat.use, c(covs, fe.use))

#Estimate
fit.tab.c5.4 = glm.nb(eq, df.use, control = glm.control(maxit = 20))
se.tab.c5.4 = fit.tab.c5.4 %>%
  coeftest(.,
           vcov = vcovHC(., type = "HC0", cluster = df.use$iso2c))

#-----------------------
#(7) Combine
#-----------------------
stargazer(list(fit.tab.c5.1,
               fit.tab.c5.2,
               fit.tab.c5.3,
               fit.tab.c5.4),
          se = list(se.tab.c5.1[,2],
                    se.tab.c5.2[,2],
                    se.tab.c5.3[,2],
                    se.tab.c5.4[,2]),
          omit = c("iso2c", "year", "Constant"),
          type = "text")
